An Algorithm for the Matrix Lambert W Function

An algorithm is proposed for computing primary matrix Lambert W functions of a square matrix A, which are solutions of the matrix equation WeW = A. The algorithm employs the Schur decomposition and blocks the triangular form in such a way that Newton’s method can be used on each diagonal block, with a starting matrix depending on the block. A natural simplification of Newton’s method for the Lambert W function is shown to be numerically unstable. By reorganizing the iteration a new Newton variant is constructed that is proved to be numerically stable. Numerical experiments demonstrate that the algorithm is able to compute the branches of the matrix Lambert W function in a numerically reliable way.